A coarse grained model for granular compaction and relaxation 
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We introduce a theoretical model for the compaction of granular materials by discrete vibrations 
which is expected to hold when the intensity of vibration is low. The dynamical unit is taken to 
be clusters of granules that belong to the same collective structure. We rigourously construct the 
model from first principles and show that numerical solutions compare favourably with a range of 
experimental results. This includes the logarithmic relaxation towards a statistical steady state, the 
effect of varying the intensity of vibration resulting in a so-called "annealing" curve, and the power 
spectrum of density fluctuations in the steady state itself. A mean field version of the model is 
introduced which shares many features with the exact model and is open to quantitative analysis. 

PACS numbers: 05.40. +j. 46.10.+Z, 64.60.L, 81.05. Rm, 81.20.Ev 



I. INTRODUCTION 



Extrapolating bulk properties from the underlying mi- 
croscopic dynamics is generally more difficult with gran- 
ular materials than with g difficulty that has been 
attributed, at least in part, to the lack of thermal av- 
eraging Unlike molecules, granules are static at 
room temperature and so cannot explore phase space 
without some external impetus. For example, consider 
a column of loosely packed granules in a cylindrical con- 
tainer, where loosely packed means that there are typi- 
cally large gaps or voids between neighbouring granules. 
It is energetically favourable for the granules to collec- 
tively reorganise to a state which minimises these voids, 
since a more compact column will have a lower centre 
of gravity and hence a lower potential energy. That this 
does not occur spontaneously is a direct consequence of 
the lack of thermal motion. One way to allow the column 
to evolve is simply to tap or otherwise perturb the con- 
tainer, thus giving the granules a small amount of kinetic 
energy with which to rearrange. This process has been 
studied empirically in the context of industrial applica- 
tions 1^], but only recently have attempts been made to 
try to understand the fundamental dynamics involved. 

Mehta et.al. [Q-^ employed a non-sequential Monte 
Carlo algorithm to simulate the process on a microscopic 
level. Non-sequential means that granules are allowed to 
move and settle simultaneously, which is important in 
this context since it allows for the cooperative reorgan- 
isation of granule-granule contacts. These simulations 
predict that granular media should relax on two time 
scales, corresponding to individual granule motion and 
collective processes respectively. However, this is not in 
accord with the experimental work of Knight et.al. 0. 
They measured the rate of compaction in a column of 
monodisperse glass beads that was subjected to discrete 
vertical vibrations. The plot of density against the num- 
ber of vibrations was found to be best described by 
p{t) ~ {\ogt)~^, where the time ordinate t is proportional 
to the number of taps. One possible reason for the dis- 



crepancy between the simulations and the experiments 
may simply be that the regimes of vibration intensity 
studied were different. The smallest vibration considered 
in the simulations corresponds to a 5% increase in volume 
at every tap, which is much more than the experiments 
involved. 

A number of models embracing a variety of theoreti- 
cal approaches have been introduced to try and account 
for the experimental findings. Of those we are aware of, 
one is a phenomenological macroscopic model but 
the remainder are all microscopic in nature. The slow 
relaxation has been attributed by Ben-Naim et.al. to 
the large number of reconfigurations required to bring 
enough small voids together to make one void large 
enough to absorb another granule . de Gennes also 

chose to focus on the voids and found that a Poisson 
distribution of void sizes could give rise to the expected 
inverse logarithmic relaxation Coglioti et.al. have 

introduced a lattice model in which each granule can be 
in one of two states with each state corresponding to a 
different geometrical orientation ||l^,|l^ . The motion be- 
tween neighbouring granules is constrained by their rel- 
ative orientations, hence the rate of relaxation in their 
model is governed by a form of geometrical frustration. 

In this paper, we introduce a model for granular com- 
paction which is neither macroscopic nor microscopic but 
instead lies somewhere between these two extremes. It 
is coarse grained in that it takes clusters of granules as 
its dynamical unit rather than individual granules. This 
approach is based on the picture of granular interactions 
described by Mehta et.al. in relation to their simula- 
tions 1^-^, except that here we are interested in the 
limit of weak vibrations. The resulting model is strik- 
ingly similar to one already devised by Bak and Snep- 
pen in a wildly different context, that of biological evo- 
lution 16|. In Sec. |l| the model is described in detail 
and its physical basis is explained. Careful considera- 
tion is given to the range of validity of our assumptions. 
Results of numerical simulations are compared to the ex- 
perimental findings in Sec. |III| . The exact solution of a 
mean field version of the model is investigated in Sec. IV. 
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Finally, with give a summary of the model in Sec. 



II. THE MODEL 

Mehta et.al. picture the granular media as being sub- 
divided into local clusters, as in Fig. |l|(a), where a cluster 
is defined as a group of granules belonging to the same 
multi-particle potential well A vibration with an 

intensity equivalent to the binding energy of a granule to 
its well causes that granule to be ejected and move in- 
dependently of the others. Under weaker vibrations, all 
the granules remain in the well but still reorganise collec- 
tively, albeit on a slower time scale to individual particle 
motion. Although this description seems to be valid for 
the range of intensities of vibration considered in their 
simulations, it clearly fails for the much lower intensi- 
ties relevant to the experiments We believe that 
the picture is essentially correct but needs to be modi- 
fied to describe the behaviour of the system deep in the 
collective relaxation regime. To do this, we first need to 
closely analyse exactly what is meant by a multi-particle 
potential well. 

Any given configuration of an ensemble of particles can 
be represented by a single point in the space of all possi- 
ble configurations. Each allowed configuration has a well 
defined potential energy, and so the time evolution of the 
ensemble under gravity can be described by a walk in con- 
figuration space over a potential energy landscape. Now, 
the preferred state for each individual granule is simply 
resting at the bottom of the container. If the granules did 
not interact, then the ensemble would trivially evolve to 
the global minimum with every granule in its preferred 
state, i.e. all resting on the bottom. Of course, real 
granules do interact, and one granule moving downwards 
will inevitably push some of the surrounding granules 
upwards slightly. The ensemble is thus frustrated in that 
it cannot simultaneously satisfy each granule's tendency 
to move downwards. In terms of the potential energy 
landscape, this frustration results in a rugged landscape 
with many local minima separated by barriers of various 
heights. A schematic example is given in Fig. ||, where 
for clarity we have compressed the entire configuration 
space onto a single axis. 

The ensemble will be at a local minimum between per- 
turbations. The effect of the perturbation is to move 
the ensemble to a point higher up on the landscape be- 
fore it again relaxes, possibly to a different minimum. 
For the low-energy perturbations we are concerned with 
here, the ensemble will usually move between nearby 
minima and consequently only a small number of gran- 
ules will change their position or orientation. Following 
Mehta et. al. we assume that these granules typically be- 
long to some sort of collective structure, such as an arch 
or bridge. Thus the system can be subdivided into lo- 
calised clusters, where a cluster is now defined as the 
unit of collective reconfiguration. Furthermore, we map 



the system onto a lattice in which every site corresponds 
to a single cluster, as in Fig. |l|. This lattice representa- 
tion is implicitly static and so will not be valid if there is 
any form of global motion in the system, such as convec- 
tion or surface flow, although it should still hold if there 
is only a limited amount of local motion. Large per- 
turbations will involve reorganisation on a system-wide 
scale and the rapid rearrangement of cluster boundaries, 
so the lattice representation is again expected to fail in 
such situations. 
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FIG. 1. An example of the process of subdividing granu- 
lar media into local clusters, given here for the case of two 
dimensions, (a) A collection of circular granules separated 
into clusters. The thick lines represent boundaries between 
neighbouring clusters, (b) The corresponding lattice repre- 
sentation. Each site denotes a single cluster. 

We now have a lattice of clusters, each of which move 
on their own individual potential energy landscapes. 
During the perturbation, each cluster is kicked to a point 
higher up on its landscape, and those that subsequently 
relax to a new minimum have collectively reconfigured. 
When a cluster reconfigures the contacts between it and 
adjacent clusters will be redistributed in a highly non- 
trivial manner, the pattern of stress lines will be locally 
distorted and the boundaries between adjacent clusters 
may shift slightly to accommodate different granules. As 
a consequence, there will be a significant change in the 
landscapes of the cluster itself and those near to it. In 
particular, we note that the heights of barriers between 
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minima will change. It may seem possible for one of the 
nearby clusters to move a significant distance on its new 
landscape before finding a minimum, effectively consti- 
tuting another reconfiguration event. However, this con- 
tradicts the definition of a cluster as the fundamental unit 
of collective reconfiguration, since any two clusters that 
interact in this way should have been treated as a single 
cluster in the first place. Thus it can safely be assumed 
that nearby clusters will not reconfigure, although the 
heights of barriers in their landscapes will still change. 

Significant progress can be made if we do away with 
the landscapes altogether and just deal with the heights 
of barriers between minima instead. Indeed, as we are 
only interested in the limit of weak perturbations, we 
can go one step further and disregard all but the small- 
est barrier, since this will almost always be the one that is 
involved anyway. Each reconfiguration is assumed to al- 
ter the landscapes in such a complicated manner that, to 
good approximation, the height of a barrier can be taken 
to be a random number drawn from a suitable probabil- 
ity distribution. Although this distribution is in general 
unknowable, we have found the model to be robust to a 
variety of different choices, including uniform, exponen- 
tial and Gaussian (robustness means that the essential 
behaviour of the system remains unchanged with respect 
to the modifications tried). We subsequently use the uni- 
form probability distribution P{b) for barrier height 6, 
where 



P(b) = <! ^ ^ S [O'l]' 
^ ' I otherwise. 



(1) 



Consider now the effect of the external perturbation on 
just a single cluster with a barrier height of 6ciust- Sup- 
pose that the effect of the perturbation is for the cluster 
to gain an energy of er and to move to a corresponding 
point higher up on its landscape. If er < ^ciust, the clus- 
ter cannot cross even its lowest barrier and so we can be 
sure that it will relax to the same minimum that it was 
at before. However, if er > 6ciust then there is a non-zero 
probability that the cluster will reconfigure. We take this 
probability to be of the form 




-^1 for er > feciust, 
for er < fociust, 



(2) 



where /i is a dimensionless constant. This may appear to 
be a somewhat arbitrary choice, but a number of varia- 
tions with a suitable cut-off at er = &ciust were tried, and 
no essential difference in system behaviour was observed. 
The choice of (|^) was made since it is exponential in form, 
implying some sort of underlying Poisson process, and it 
has the correct asymptotics for er ^ciust and er — * cxd. 




current configuration 



Configuration 

FIG. 2. Schematic example of a potential energy landscape 
for an ensemble of granules in configuration space. The en- 
semble currently lies at the local minimum marked. The 
smallest barrier to an adjacent minimum has a height of feciust- 

When the container is vibrated, the associated energy 
impulse is distributed in some undefined manner to all 
the clusters in the system. We have observed little qual- 
itative difference arising from distributing this energy 
stochastically and henceforth assume that each cluster 
receives the same energy er- It should be clear from (j^) 
that the cluster with the smallest barrier in the system, 
say of height 6min, is the most likely to reconfigure. With 
this observation, we can make a further simplification 
that also makes little difference to the system behaviour, 
which is to assume that the cluster that reconfigures first 
is always the one with the barrier height of 5min. Thus 
there is no longer any need to simulate every perturbation 
until the cluster reconfigures, we can instead just recon- 
figure the cluster immediately and advance the time by 
an amount St, where 



St oc exp < /i 



cr 



er 



(3) 



This is the expected number of perturbations of energy 
er required until the cluster with barrier height femin re- 
configures, and is the reciprocal of (||). For 6niin < er, St 
is taken to be infinite. 

We are now in a position to describe the model al- 
gorithmically. The granular media is represented by a 
lattice, each site of which corresponds to a unit of collec- 
tive reconfiguration, ie. a cluster. The model is robust to 
variations in lattice connectivity, so without loss of gener- 
ality we choose a simple cubic array. Each cluster {i, j, k) 
has an associated potential energy barrier against recon- 
figuration, bijk, drawn from the probability distribution 
P{b) given in (|l|). The external perturbation takes the 
form of an energy impulse distributed uniformly through- 
out the system, each cluster receiving an amount er. At 
each algorithm step, the cluster with the smallest barrier 
in the system, fomin, is found. If er < fomin then the per- 
turbation is too weak to cause any reconfiguration events, 
the system is frozen and the simulation is complete. If 
er > bnun, the cluster in question and the 6 clusters ad- 
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jacent to it are reconfigured - that is, their barriers are 
redrawn from the same probabiUty distribution as before. 
The real time is increased by an amount St defined in (||), 
and the simulation moves on to the next algorithm step. 
Note that wc do not employ periodic boundary condi- 
tions, instead clusters at the faces, edges or corners of 
the lattice simply have 5, 4 or 3 adjacent clusters, re- 
spectively. 

Numerical solutions of the model are presented in the 
following section. For now, we would like to remark 
upon the strong similarity between this model and a 
model of biological evolution already devised by Bak and 
Sneppen |^5|. The lattice sites in their model repre- 
sent different species, each of which is assigned a bar- 
rier against mutation corresponding to the smallest bar- 
rier between local optima on a rugged fitness landscape. 
The species mutate and interact with adjacent species in 
much the same way that clusters reconfigure and inter- 
act with adjacent clusters in our model. The primary 
difference between the models is that, whereas clusters 
cannot move higher than er on their potential energy 
landscapes, corresponding to the strength of the external 
impulse, species are subject to no such energetic con- 
straints (there is no such thing as the "conservation of 
fitness" ) and move around their fitness landscapes spon- 
taneously. As long as this difference is borne in mind, we 
can draw upon the plethora of results already accumu- 
lated for the evolution model in analysing our model of 
compaction (for a review, see |16 ). 



III. COMPARISON TO EXPERIMENTS 



We begin by describing the numerical solution of the 
model for a system comprising of N clusters. The dis- 
tribution of barrier heights, Q{b), is defined such that a 
proportion Q{b)db of the clusters have a barrier height in 
the range b to b + 5b. As the system evolves, Q{b) ex- 
hibits two qualitatively different regions, one for large b 
and one for small b. Large barriers have either not been 
touched since the simulation began, or (more likely) they 
have been redrawn from the uniform distribution P{b) 
as the consequence of an adjacent cluster reconfiguring. 
As such, Q{b) for large b must also be uniform, except 
for statistical fluctuations. The situation is more com- 
plicated for small barriers since there is now the added 
possibility of being selected as the minimum of the sys- 
tem. Very small barriers are unlikely to last long and so 
Q{b) tails off to zero as 6 ^ 0. The boundary between 
these two regions is given by the gap function G{t), which 
is the largest barrier height that has ever been the min- 
imum of the system. Finding the minimum barrier and 
giving it a new value can be viewed as a flux from the 
region b < G{t) to the region b > G{t). When there are 
no barriers left in the region b < G(t), larger barriers will 
be selected as the minimum and so G{t) will increase. 
If there were no interactions, there would only be this 



unidirectional flux and G{t) would slowly approach 1 as 
t —^ oo. However, with interactions there is also a flux 
in the reverse direction, from b > G{t) to b < G{t), cor- 
responding to the new values given to the barriers of 
adjacent clusters. Hence G{t) in fact approaches a con- 
stant value b* G (0, 1), where b* is a function of the lattice 
connectivity and the system size N. 

We have not yet considered the effect of the parame- 
ter er- This appears in the equation for St, the time step 
between successive reconfiguration events, which also de- 
pends on the current value of the minimum barrier (|^). 
It can be seen from (|^) that St becomes singular when the 
minimum barrier is greater than or equal to er . If er > b* 
then this can never happen, since the minimum fluctuates 
between and G{t), and G{t) — > 6* as t — > oo. Accord- 
ingly the system approaches a statistical steady state in 
which St fluctuates around some constant value. By con- 
trast, if er < b* then it now becomes possible for G(t), 
and hence also the minimum, to take values close to er. 
As it does so, St will diverge and the system will freeze 
into a state in which every cluster has a barrier greater 
than er and there can be no further reconfigurations. An 
example of how G(t) depends on er is given in Fig. ^ for 
a 40 X 40 X 40 lattice, for which 6* « 0.21. 
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FIG. 3. Plot of the gap function G(t) for various values of 
er, for a 40 x 40 x 40 lattice. Key: Plus signs, er = 0.4. As- 
terixes, er ~ 0.3. Open circles, er = 0.25. Crosses, er ~ 0.2. 
Filled circles, er ~ 0.15. Dots, er ~ 0.1. Note that in this 
and all subsequent plots we have taken the time step to be 
St = exp{er/(er — bmin)}/N, where A'' is the system size, so 
the units on the time axis are arbitrary. 

The model has so far been described in terms of the en- 
ergy impulse per cluster er and the barrier distribution 
Q(b). However, the experimental results were given in 
terms of an acceleration parameter F and the density p. 
Before comparing the model with the experimental re- 
sults, we must first consider how these two sets of quan- 
tities are related. We start with er and F. The accel- 
eration parameter F is defined as the peak acceleration 
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during the perturbation scaled by gravity, F = amax/ ■ 
This was also found to be the relevant parameter for the 
stability of a bead heap under vibration Although 
it seems reasonable that a higher T should mean more 
energy is distributed throughout the system and hence 
a higher er, the precise relationship is likely to be very 
complex and we have been unable to derive a formula 
relating the two. Instead we simply assume that, for the 
small vibrations considered here, the relationship is ap- 
proximately linear, er oc F. 

Trying to quantify the relationship between the barrier 
distribution and density is more problematic since a po- 
tential energy barrier is an intrinsically abstract concept. 
Nonetheless, a rough formula can be derived as follows. 
Consider an individual cluster with a barrier 6ciust and 
density Pciust- The cluster's horizontal cross sectional 
area is assumed to remain roughly constant throughout 
the compaction process, so the typical vertical separation 
between the granule centres will be inversely proportional 
to Pciust- The cluster cannot reconfigure unless this ver- 
tical separation is increased to the order of the granule 
diameter, thus allowing the granules to move over one an- 
other. Since the granule diameter is constant, the change 
in height required for reconfiguration will also depend 
inversely upon Pciust- The potential energy gained by a 
particle is, of course, proportional to its height increase, 
so 6ciust also varies inversely with Pciust- Extrapolating 
this result over the entire system amounts to finding the 
mean barrier height 6, so finally we have 

b^p-'. (4) 

This derivation is simplified in that, for instance, it does 
not incorporate the effect of adjacent clusters on the value 
of &ciust- We expect it to work for overall trends in density 
but not for small fluctuations. 

We are now in a position to test the model against the 
experimental results. As mentioned in the introduction, 
the density was experimentally found to relax inverse log- 
arithmically with time, p{t) ~ {\ogt)^^ 0. From (^) the 
corresponding relationship in terms of the mean barrier 
height is therefore b{t) ~ logi, which will show up as a 
straight line on a graph of b{t) vs \ogt. Such a graph is 
given in Fig. ^ for a range of values of er- Linear be- 
haviour is apparent over a broad range of densities for 
er > b* , confirming logarithmic relaxation towards the 
statistical steady state. For er < b*, the relaxation is 
initially logarithmic but slows down as the frozen steady 
state is approached. Note that although the logarithmic 
behaviour is robust, the actual values on the axes depend 
upon which of the various arbitrary choices mentioned in 
the previous section have been made and hence have no 
physical meaning. 

Little has been said so far about initial conditions. 
Before the first selection of the minimum barrier Q{b) 
is uniform over the entire range [0,1], so that even a 
small er will cause a significant amount of reconfigura- 
tion. This corresponds to a state of minimum compactiv- 
ity which is very difficult to attain experimentally. For 



instance, there will always be a certain amount of back- 
ground noise, and the granules added later to the appa- 
ratus will impact upon those already present, inevitably 
causing some compaction. Instead, the experiments al- 
ways started from a slightly compacted state with a den- 
sity fraction of 0.577±0.005. This initial compaction can 
be incorporated into the model by shifting the time axis 
so that the origin corresponds to when G{t) first becomes 
greater than a parameter 6init > 0. Values of er ~ foinit 
or less are too small to cause any significant further com- 
paction. This is readily apparent in Fig. ^, where we 
have plotted b in the limit t oo against er. The line 
is flat for er < ^init, increases linearly for 6init < er < 6* 
and levels out again for higher er. This should be com- 
pared with the corresponding experimental plot, which is 
Fig. 3 in , from which we estimate that b* corresponds 
to F w 3. 
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FIG. 4. b(t) vs. Int for a range of values of er. The data 
was taken from single runs on a 40 x 40 x 40 lattice, for which 
b* ~ 0.21. From top to bottom, the values of ep are: 0.4, 0.3, 
0.25, 0.2, 0.15, 0.1. Solid lines have been used for ep > b* and 
dashed lines have been used for ep < fo*. 

An apparently anomalous feature of Fig. ^ is that the 
highest densities are to be found, not for large er, as 
might be expected, but instead for values of er near the 
threshold value b* . This occurs because of finite size ef- 
fects. Recall that, for er > 6*, the barrier distribution 
evolves to a state which is uniform for b > b* with a tail 
for b < b* . It is the very existence of this tail, which 
disappears in the thermodynamic limit — > oo, that re- 
duces the mean barrier b for finite systems. When er 
is slightly less than b* then, although the uniform re- 
gion is slightly broader, the selection process can remove 
some of the barriers from the tail permanently and so 
the net effect is to increase b. An even greater degree 
of compaction can be obtained if a system with er > b* 
is first allowed to self-organise to the statistical steady 
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state, then er is slowly reduced to zero to remove as 
much of the tail as possible. Quickly reducing ep will not 
give enough time for the selection process to work be- 
fore the system froze and so b would barely change. An 
example of this process is given in Fig. ^, where to ac- 
centuate the finite size effects a 4 x 4 x 4 lattice was used. 
Nowak et.al. have produced similar plots from their ex- 
periments, which they regard as a type of annealing pro- 
cess [ pl]p^ . They label the lower branch of the graph, 
when the intensity of vibration is increased for the first 
time, as "irreversible". In the language of our model, 
we prefer to call this the self- organising branch. The 
self-organising branch meets an upper reversible branch 
around the point T* « 3. This is to be expected since, 
as mentioned in the previous paragraph, this value of T 
corresponds to the threshold value b* , that is, the point 
at which the system can self-organise into the statisti- 
cal steady state. According to the model, the change in 
density along the upper branch is due to the effects of 
finite size, so there should be a greater variation when 
larger beads are used in the same sized apparatus. This 
is in agreement with the experiments except for when 
the largest bead size was used 



18[ . In this case. 



al- 
though the overall density variation was the greatest, a 
disproportionately large amount of it occurred along the 
self-organising branch, possibly due to the cylinder walls 
aligning the beads into a highly compact crystalline con- 
figuration. Another feature observed in the experiments 
is that the threshold value F* appears to increase when F 
is updated more rapidly. The model agrees with this and 
attributes it to the larger number of steps that will take 
place before the system has had time to self-organise. 




Vibration intensity, e^^ 

FIG. 5. The mean barrier height b in the final steady state 
as a function of ep. Note that since b oc (po — P)~^ ^ P the 
vertical axis can also be identified as the (approximate) den- 
sity. The simulations were performed on a 10 x 10 x 10 lattice 
and averaged over 1000 runs, fcinit = 0.08 and b* « 0.25. 



For er > b* the steady state is statistical in nature, so 
another test for the model would be to compare the fluc- 
tuations of b around its steady state value to the fluctu- 
ations in density measured experimentally. However, as 



previously mentioned, the argument relating 6 to p is not 
expected to hold for small changes. The change in density 
caused by, say, a single reconfiguration event will be sen- 
sitive to the exact positions of a large number of granules 
at that instant in time. The experimental plot of density 
fluctuations is Gaussian in form |11|, indicative of the 
large number of independent factors involved. A more 
revealing distribution is the power spectrum of density 
fluctuations, S{f), where the frequency / is measured 
in units of (taps)~^. Experimentally, S{f) was found to 
obey the power law S{f) oc f'^ , with 5 = 0.9 ± 0.2, for 
a broad range of /. Apart from finite size effects, the 
model predicts a power law with 5 = 1 . When large 
intensities of vibration were applied in the experiments, 
the power law behaviour was broken up by regions with 
(5 = 0, 0.5 or 2. We cannot account for this and attribute 
it to the expected breakdown of the model for large vi- 
brations. 
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Vibration Intensity, e^, 

FIG. 6. Annealing curve for a 4 x 4 x 4 lattice, for which 
b* ~ 0.38. er first increases from 0.04 to 0.68 in steps of 
0.04 (filled circles), then decreases by the same step size from 
0.68 to 0.04 (open circles). Finally, er is increased up to 0.68 
again (asterixes). &init was set at 0.08. Each simulation was 
run until t ~ 156, and the final plot was averaged over 1000 
such runs. 



We end this section by briefly considering how the 
model might also be applied to a set of related experi- 
ments. Jaeger et.al. ||l^ have shown that the angle of 
repose 6(t) of granular media in a half- filled cylindri- 
cal drum relaxes according to 9{t) ^ logi when vibrated. 
Furthermore, they also demonstrated the existence of a 
threshold in the intensity of vibration below which the 
relaxation was qualitatively slower. If we ignore the 
compaction process which presumably occurs in the bulk 
of the pile, then the typical vertical separation between 
granule centres is now proportional to tanO, although for 
the range of angles involved we can use tan 9 9 instead. 
We can now repeat the argument given earlier for density 
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and find that the relationship between the mean barrier 
and the slope is b ~ 9, to first order. Hence the model 
also predicts relaxation of the form 9{t) ~ logt and the 
existence of the threshold in the intensity of vibration. 
However, we have reservations in applying the model to 
this new geometry since it blatantly involves a global, al- 
beit slow, movement of granules over the surface, some- 
thing which we have explicitly stated the model does not 
cater for. It should also be mentioned that other theo- 
retical explanations for this behaviour have already been 
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IV. MEAN-FIELD ANALYSIS 

The picture presented thus far can be extended by 
considering a mean field version of the model which is 
open to quantitative analysis. This simplified model ex- 
hibits many of the traits apparent in the exact model, 
especially in the relaxation towards the statistical steady 
state. However, it behaves very differently in the steady 
state itself, and we refer the reader elsewhere for analysis 
of the original model in this much studied regime ||l5|,|l6| . 
The required mean field approximation is to be achieved 
in two stages. First, all spatial definition is removed. 
This means that, when the cluster with the smallest bar- 
rier in a system of N clusters is found and reconfigured, 
K other clusters are chosen at random from the remain- 
ing A'^ — 1 and their barriers given new values. These 
K clusters are equivalent to the adjacent clusters in the 
original model, so for example K ~ 6 corresponds to a 3 
dimensional system. The second simplification is to as- 
sume that N is very large. In this way the system can be 
described by continuous rather than discrete variables, 
to within an error margin of 0{1/N). 

For the first part of this section, the evolution of the 
system will be described in terms of a time variable r 
which increases by 1/iV between successive reconfigura- 
tions. The inclusion of the variable time step given in (||) 
will be postponed until later. The system is described by 
the cumulative barrier distribution C(6,t), which is de- 
fined as the proportion of clusters with barriers less than 
6 at time r and is related to (5(6, r) by 



C(6,t) = / Q{x,T)dx . 
Jo 



(5) 



The time scale has been normalised to one reconfigura- 
tion per cluster per unit r, so C(&, r) evolves according 
to 



dC{b,T) 
dr 



-0(6 - 6„in(r)) - K C{b, r) + 6 (X + 1) , 



(6) 



where 6min(''') is the value of the minimum barrier in 
the system at time r and 9{b) = 1 for 6 > and oth- 
erwise. The removal of the minimum barrier has the 



effect of reducing C(6, r) for all values of 6 > 6niin('r) but 
leaves it unchanged for 6 < 6min (''')■ This is handled by 
the first term on the right hand side of (||). In a simi- 
lar manner, the second and third terms account for the 
selection of the K random nearest neighbours and the 
K + 1 new barrier values, respectively. It is straightfor- 
ward to check that (||) preserves C(0, t) = 0, C(l, r) = 1 
and C(6i,t) > C(62,t) for 6i > 62, for all values of r. 

The rate equation (g) is not yet in a closed form be- 
cause it involves the unknown quantity 6min(T). We 
might naively try to write down a second equation giv- 
ing 6ijiin(r) in terms of C(6, t), perhaps something like 
C(6niin('''), ''■) — l/N. However, it must be recalled that 
errors of 0{1/N) have already been made in going from 
the discrete model to this continuous description, and so 
C (6, r) cannot be used to this degree of accuracy. In- 
deed, any attempt to define the minimum barrier within 
a continuum framework is doomed to failure for this very 
reason. We are forced to conclude that there can be no 
set of closed equations in terms of C(6, t). All is not lost, 
however, since this problem can be partially circumnavi- 
gated by use of the gap function G{t). As before, G{t) is 
defined as the highest value that 6miii(T) has ever taken, 
or more formally. 



G(r) 



sup b„ 

0<Z<T 



(7) 



Values of 6 greater than G(t) must by definition be 
greater than every value 6min has taken up to a time 
T. This allows for (g) to be simplified to 



dCjb, t) 
dr 



-{KC{b,T) + l) + b{K + l) 



(8) 



for 6 > G(t). This can be solved by substituting 
C{b,T) = a{T)b + (3{t) and comparing coefficients of 6. 
With the initial condition (7(6,0) = 6 (so 6init = 0), the 
result is 



C(6,r) 



1 



K 



(1 



-Kt 



(9) 



The fact that (7(6, t) is linear means that the barrier dis- 
tribution (3(6, t) is uniform for 6 > G(t), as expected. 
The solution (g) holds from 6=1 down to 6 « G(r) , 
where the continuum approximation starts to break down 
and we have entered into the asymptotic tail. Since there 
are only 0(l/iV) clusters in this tail, the value of G(r) 
will correspond to the point at which G(6, r) is zero, 
ie. C{G{t),t) — 0. Together with (^) this allows for the 
time dependent form of G(r) to be found. 



G(r) 



1 



-Kt 



K + l-e 



-Kt 



(10) 



Ray and Jan have also found this result by an alternative 
method ||2^. The threshold value of 6 in this mean field 
model is therefore 
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b* = lim G(t) 



K 



(11) 



which is smaller than in the exact model. 

In this approximation, the mean barrier height b be- 
haves in the same way as the gap function. This 
is because, to 0{1/N), there is no tail for b < G{t) 
and the barrier distribution is uniform for b > G(r), 
so 6(t) = (1 + G{t))/2, which is just a linear rescaling. 
Hence we expect G(t) to vary logarithmically with r. 
When the expression for G{t) given in (0) is plotted 
against log r it exhibits a linear region similar to the ex- 
act model, but not extending quite as close to the steady 
state. The gradient of G{t) in this log-linear plot is 



dG'(T) dG(r) 

— T 

d(lnT) dr 



rG'(r) 



(12) 



The linear region occurs around the point where the gra- 
dient is stationary, ie. when the second derivative is zero. 



d / dG( 
d(lnr) \d{l: 



^) =r(G'(r)+rG"(r)) = 0. 



(13) 



The solution with t — corresponds to the singularity 
in Inr and can be ignored. Using (|l^), the non-trivial 
solution is 



T= — tanhy(T-f To) , 



(14) 



where the constant tq = {\n{K + 1))/K. Since the slope 
is roughly constant in this region there is no need to find 
the exact value of t that satisfies (^4|). Instead we ob- 
serve that, for large K, the tanh function is roughly equal 
to 1 for all r > 0, so an approximate solution is t « 1/K 
and hence the slope is 



dG(T) 



d(lnr) 



Ke 



(15) 



We now turn to consider the effect of the variable 
timestep St as defined in (|3|), which depends on femin 
and er- The quantity bmhi is unknown, but we know from 
the discrete model that it fluctuates between and G(t) 
and therefore substituting G(t) for 6min(''") gives a quali- 
tatively identical solution. The new time scale is denoted 
by ^(t) and is defined by 



dt 

— = exp 
dr 



er 



er - G(r) 



(16) 



For small r, G(r) = t + 0{t^) and (|l|) can be solved 
with the initial condition t(0) =0 to give 



t{r) 



(17) 



which is linear up to r = O(ep^). The behaviour of tij) 
for large r depends upon whether er is greater than. 



less than or equal to the threshold value b* = -ff^rj-. 
er > 6*, G(r) T^+i as r ^ oo and consequently 



t ~ r exp < /i 




For 



(18) 



The time scale is stretched by a constant factor, but oth- 
erwise the system approaches the same statistical steady 
state as before. For er < 6*, (16) becomes singular at the 
point T = Tcrit at which Giji^^it) — er. Since 5t diverges 
there are no more reconfigurations and the system is in a 
frozen steady state. The precise nature of this singularity 
can be found by substituting r — Td-it — e into (p^), with 
e small and positive. As e ^ 0, i(r) diverges according 
to 



dt 
d7 



where the constant 
A = /i 



er 



(19) 



(20) 



(l-er)(l-(i^ + l)er) ' 
Finally, for er = b* (n6h can be algebraically reduced to 



di 

dr 



exp < /i 



K + 1 



K 



(21) 



for large r, which is divergent. 

Now that we have confirmed that er has the same ef- 
fect in the mean field model as in the exact model, we 
need to see what it does to the rate of logarithmic decay. 
This is straightforward for er 3> b* since 



t = e^'T + 0{e^^) , 



(22) 



so to first order in e^^ the time scale is just stretched 
by a constant factor, which does not alter the gradient in 
a log-linear plot. This means that slope of G{t) vs \ogt 
is the same as the slope of G(r) vs logr and ( [Tsl ) can 
be used without modification. For instance, in the exact 
system with large er the slope is approximately 0.048 
in 3 dimensions, whereas the value predicted by ( |l5|) for 
= 6 is 0.050. 

Modifying (|l5|) to incorporate er < oo is troublesome 
and we have been unable to derive a general formula. 
Nonetheless there is still some hint of a correspon- 
dence between this analysis and the experiments. In 0| 
Knight et. al. introduce a parameter r which we call roxp 
so as not to confuse it with our r. rcxp gives a rough 
measure of the time scale of the relaxation process. We 
tentatively equate this to the quantity di/dr, and in- 
deed the experimental plot of roxp vs F looks similar to 



the form of dt/dr given in (p_6[). However, this is not 
a robust feature of the model and so it is impossible to 
come to any concrete conclusions. The experimental data 
also shows a noticeable change in behaviour for small F. 
This could be caused the system entering into the frozen 
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steady state before the logarithmic relaxation has had 
a chance to take hold, ie. when Tdit ^ although it 
could just be the effect of the initial compaction. 

Finally, we demonstrate how this analysis can be ex- 
tended to incorporate energy dissipated by a reconfigur- 
ing cluster to its nearest neighbours. Suppose that each 
adjacent cluster receives an energy ediss and immediately 
reconfigures if its barrier is smaller than this, dissipating 
a further energy Cdiss to each of its neighbours, and so on. 
Using the same mean field approximations as before, the 
net effect of this avalanche between perturbations is to 
increase the number of barriers that change value at each 
time step. Of the K random nearest neighbours, iiTediss 
will immediately reconfigure and so the total number of 
new barriers per time step dr is now 



K + ediss) + /^(i^ediss)' + K{Kedissf 



K 



1 - A'cdii 



(23) 



for ediss < • Larger values of ediss are unphysical since 
they result in perpetual reconfiguration. The new rate 
equation for C(6, r) is 



dC{h, t) 
dr 



= -9{b - 6„,in(r)) - 



K 



1 



K 



1-Ked 

which can be solved as before to give 



(24) 



C(6,t) =5+^(f-ifedi. 



1 — exp 



Kt 



1 - Ke^,, 



(25) 



for b > G{t). This is the same as the solution al- 
ready given in (^ except that K has been replaced 
by the effective number of random nearest neighbours 
K/{1 — K ediss)- The time scale is similarly stretched by 



the constant factor 1 — KeA 



Hence the inclusion of 



energy dissipation in this manner does not alter the be- 
haviour of the system, nor does it change the slope of 
G{t) in a log-linear graph. 



V. SUMMARY AND DISCUSSION 

We have presented a theoretical model for the com- 
paction of granular materials by low intensity perturba- 
tions which appears to agree well with a range of experi- 
mental results. This includes the logarithmic relaxation, 
the effect of varying the intensity of vibration resulting 
in a so-called "annealing" curve, and the power spec- 
trum of density fluctuations in the steady state. We have 
segmented the granular media into local subsystems or 



clusters which represent ensembles of granules that col- 
lectively reconfigure. Associated with each cluster is a 
potential energy barrier against reconfiguration. When- 
ever a perturbation gives a cluster enough energy to cross 
over its barrier into a new configuration, nearby clusters 
are disrupted and their barriers take on new values. The 
system behaviour is dominated by this dynamical inter- 
action between clusters and fine detail such as the choice 
of distribution for the barrier values makes little or no 
difference. Indeed, it is this very robustness that leads 
us to hope that the model might correctly describe the 
mechanism underlying the compaction process, despite 
its algorithmic simplicity. 



It has been suggested that standard statistical mechan- 
ics can be applied to granular materials if the fundamen- 
tal quantities involved are suitably reinterpreted pl| , p2[ . 
Volume plays the role of energy, and the quantity con- 
jugate to volume is known as compactivity, which is the 
analogue of temperature. The compactivity is infinite 
when the system is at its maximum volume and zero 
when it is at its minimum. Our model can also be de- 
scribed in terms of volume rather than energy since the 
external perturbations increase the volume of the sys- 
tem as well as its energy. Hence we can assign a volume 
barrier to each cluster which must be exceeded for re- 
configuration to take place. In this way, we can see the 
beginnings of a link to the modified statistical mechanics, 
perhaps with the barriers being in some way related to 
the compactivity. This is just speculation, however, and 
further investigation is required. There are also be many 
ways in which the model can be enhanced make it more 
physically realistic. For instance, the model is currently 
isotropic, but real granular media exhibits a density gra- 
dient with the densest regions near the bottom. 



There is another way to compact granules into a 
smaller volume, and that is simply to apply a uniform 
pressure. This forces the granules to rearrange into a 
higher density state, as with the perturbation-induced 
compaction studied in this paper, although the granules 
are now also subject to deformation and fracturing. A 
theoretical model for compaction by applied pressure has 
been proposed which treats the media as being comprised 
of a number of subsystems, each of which is associated 
with a pressure barrier |23). This obviously bears some 
similarity to the approach we have adopted in construct- 
ing our model. A crucial difference is that the subsystems 
in the pressure model do not interact and the values for 
the barriers are simply drawn from a suitable distribu- 
tion. In our model, the choice of distribution is unimpor- 
tant and it is the dynamical interactions between sub- 
systems that dominates the system behaviour. It would 
be interesting to see if the interacting cluster picture can 
be applied to this or any other experimental situation 
involving granular materials. 
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